{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7f3912d7-38a3-491a-a2dd-4c58768437e8",
   "metadata": {},
   "source": [
    "<span style=\"font-size:20px;\">【开源实习】QAOA量子优势的衡量</span>\n",
    "\n",
    "<span style=\"font-size:17px;\">1 任务背景</span>\n",
    "\n",
    "QAOA是常用的量子近似优化算法，它可以看作是在经典平均场算法上添加量子涨落。而为了衡量QAOA在哪些问题上能展现量子优势，可以使用《Mean-Field Approximate Optimization Algorithm》文中提出的Lyapunov指数（文中式70）进行衡量。用随机SK模型作为例子，用Lyapunov指数来衡量QAOA的量子优势，最终复现FIG.7的四幅图\n",
    "\n",
    "<span style=\"font-size:17px;\">2 任务需求及实现方案</span>\n",
    "\n",
    "基于mindquantum的Hamiltonian构造哈密顿量来求解QAOA的基态能量E0。Mean-Field AOA算法的能量和Lyapunov指数只需通过矩阵运算即可得到 最终用matplotlib输出关于时间的函数E(t)即可\n",
    "\n",
    "<span style=\"font-size:17px;\">3 算法实现及效果</span>\n",
    "\n",
    "下面基于mindquantum来实现该算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cd9d8d1c-2c28-4b1b-9141-1723ae43d42a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import expm, eigvals\n",
    "from mindquantum.core.gates import X\n",
    "import matplotlib.pyplot as plt\n",
    "from mindquantum import Hamiltonian\n",
    "from mindquantum.core.operators import QubitOperator\n",
    "from scipy.sparse.linalg import eigsh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bbe149f-898b-47c3-a0e4-fd89b86f8672",
   "metadata": {},
   "source": [
    "这是两个输出最终结果图的函数\n",
    "\n",
    "第一个函数输出的是FIG7.3和FIG7.4,横坐标是归一化的时间s，纵坐标是Lyapunov指数，最终输出了Mean-field AOA算法的前三个Lyapunov指数\n",
    "\n",
    "第二个函数输出的是FIG7.1和FIG7.2,横坐标是归一化的时间，纵坐标是E-E0，QAOA算法的20条蓝线输出的是QAOA算法前第2-21个能量E-E0 红线输出的是Mean-field AOA算法的能量E-E0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d5915105-3369-4234-b10f-89dcb1881b93",
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_graph(lamda_first, lamda_second, lamda_third):\n",
    "    \"\"\"\n",
    "    Show the top three lyapunov exponents over time.\n",
    "\n",
    "    Args:\n",
    "        lamda_first: The largest Lyapunov exponents.\n",
    "        lamda_second: The second largest Lyapunov exponents.\n",
    "        lamda_third: The third largest Lyapunov exponents.\n",
    "    \"\"\"\n",
    "    x_label = [i / len(lamda_first) for i in range(len(lamda_first))]\n",
    "    plt.figure(figsize=(8, 8))\n",
    "    plt.plot(x_label, lamda_first, label=\"first\", color=\"red\")\n",
    "    plt.plot(x_label, lamda_second, label=\"second\", color=\"gray\")\n",
    "    plt.plot(x_label, lamda_third, label=\"third\", color=\"gray\")\n",
    "    plt.xlim(0, 1)\n",
    "    plt.ylim(0, 1)\n",
    "    plt.xlabel(\"s\")\n",
    "    plt.ylabel(\"λ\")\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def draw_graph1(eig_all, E_meanfield):\n",
    "    \"\"\"\n",
    "    Show the top twenty exact eigen spectrum and the energy returned by the mean-field AOA over time.\n",
    "\n",
    "    Args:\n",
    "        eig_all: the exact eigen spectrum.\n",
    "        E_meanfield: the energy returned by the mean-field AOA.\n",
    "    \"\"\"\n",
    "    x_label = [i / eig_all.shape[0] for i in range(eig_all.shape[0])]\n",
    "    plt.figure(figsize=(8, 8))\n",
    "    for i in range(20):\n",
    "        plt.plot(x_label, eig_all[:, i + 1], color=\"blue\")\n",
    "    plt.plot(x_label, E_meanfield, color=\"red\", linestyle=\"--\")\n",
    "    plt.xlim(0, 1)\n",
    "    plt.ylim(0, 4)\n",
    "    plt.xlabel(\"s\")\n",
    "    plt.ylabel(\"E-E0\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0a6c2fbf-d5bd-4f73-ae62-10d14df62ad6",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Problem:\n",
    "    def __init__(self, num_layers, couplings, driver=X):\n",
    "        self.num_qubits = len(couplings) - 1\n",
    "        self.num_layers = num_layers\n",
    "        self.local_fields = couplings[:-1, -1]\n",
    "        self.couplings = couplings[:-1, :-1]\n",
    "        self.driver = driver"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bcc315c-c621-497e-8f7d-eefc63926809",
   "metadata": {},
   "source": [
    "这个函数实现了量子涨落矩阵的计算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bb3e31fa-4975-41d7-818d-a8d78fdcccb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fluctuation_matrix(problem, S, solution, beta, gamma):\n",
    "    \"\"\"\n",
    "    Returns the Gaussian fluctuation matrix for a given point along the mean-field trajectories.\n",
    "\n",
    "    Args：\n",
    "        problem (Problem)\n",
    "        S (np.ndarray), the initial vector of all spin vectors.\n",
    "        solution (List), the solution vector.\n",
    "        beta (List), the vector of QAOA driver parameters.\n",
    "        gamma (List), the vector of QAOA problem parameters.\n",
    "\n",
    "    Returns：\n",
    "        matrix, the fluctuation_matrix L.\n",
    "    \"\"\"\n",
    "    num_qubits = problem.num_qubits\n",
    "    local_fields = problem.local_fields\n",
    "    couplings = problem.couplings\n",
    "    A = np.zeros((num_qubits, num_qubits), dtype=complex)\n",
    "    B = np.zeros((num_qubits, num_qubits), dtype=complex)\n",
    "    t_3 = np.diag(np.concatenate([np.ones(num_qubits), -np.ones(num_qubits)]))\n",
    "\n",
    "    def n_ij_pm(idx, pm):\n",
    "        return solution[idx] * S[idx, 0] + pm * 1j * S[idx, 1]\n",
    "\n",
    "    for i in range(num_qubits):\n",
    "        A[i, i] = (\n",
    "            -beta * S[i, 0] / (1 + solution[i] * S[i, 2])\n",
    "            - gamma\n",
    "            * solution[i]\n",
    "            * magnetization(S, local_fields, couplings, num_qubits)[i]\n",
    "        )\n",
    "        for j in range(i + 1, num_qubits):\n",
    "            A[i, j] = gamma * couplings[i, j] * n_ij_pm(i, 1) * n_ij_pm(j, -1)\n",
    "            B[i, j] = gamma * couplings[i, j] * n_ij_pm(i, 1) * n_ij_pm(j, 1)\n",
    "    A += np.conj(A).T\n",
    "    B += B.T\n",
    "    top_block = np.hstack((A, B))\n",
    "    bottom_block = np.hstack((B.conj().T, A.conj()))\n",
    "    L = t_3 @ np.vstack((top_block, bottom_block))\n",
    "    return L"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dce6c07-d77c-47d0-83b0-9c792d38eb1e",
   "metadata": {},
   "source": [
    "这个函数利用上面的量子涨落矩阵来计算出lyapunov指数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f6978af5-a220-4bc4-9e30-5aad170aae1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def evolve_fluctuations(problem, t, beta, gamma):\n",
    "    \"\"\"\n",
    "    Calculate the evolution of Gaussian fluctuations over time.\n",
    "\n",
    "    Args：\n",
    "        problem(Problem)\n",
    "        t (float), the time-step of the considered annealing schedule.\n",
    "        beta (List), the vector of QAOA driver parameters.\n",
    "        gamma (List), the vector of QAOA problem parameters.\n",
    "\n",
    "    Returns：\n",
    "        The Lyapunov exponents characterizing the dynamics of the Gaussian fluctuations.\n",
    "    \"\"\"\n",
    "    num_qubits = problem.num_qubits\n",
    "    local_fields = problem.local_fields\n",
    "    couplings = problem.couplings\n",
    "    num_layers = problem.num_layers\n",
    "    Lyapunov_exponent = np.zeros((num_layers, 2 * num_qubits))\n",
    "    S = np.zeros((num_layers + 1, num_qubits, 3), dtype=np.float64)\n",
    "    S[0, :, :] = np.array([[1.0, 0.0, 0.0] for _ in range(num_qubits)])\n",
    "    S = evolve(S, local_fields, couplings, beta, gamma)\n",
    "    solutions = mean_field_solution(S[-1, :, :])\n",
    "    m = np.eye(2 * num_qubits, dtype=np.complex128)\n",
    "    for k in range(len(gamma)):\n",
    "        L = fluctuation_matrix(problem, S[k, :, :], solutions, beta[k], gamma[k])\n",
    "        m = expm(-1j * t * L) @ m\n",
    "        lyapunov_exponential_eig = eigvals(m @ m.conj().T)\n",
    "        Lyapunov_exponent[k] = np.sort(np.real(np.log(lyapunov_exponential_eig) / 2))\n",
    "    return Lyapunov_exponent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "20fcba55-374b-46fc-8ba3-65ffa1463d92",
   "metadata": {},
   "outputs": [],
   "source": [
    "def magnetization(S, h, J, num_qubits):\n",
    "    \"\"\"Returns the magnetization vector belonging to the input 'S'.\"\"\"\n",
    "    return np.array(h) + np.array(\n",
    "        [sum(J[i, j] * S[j, 2] for j in range(num_qubits)) for i in range(num_qubits)]\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3bea2f31-1705-4d14-b26b-0fd40ff64327",
   "metadata": {},
   "outputs": [],
   "source": [
    "def V_P(alpha):\n",
    "    \"\"\"Returns the rotation matrix resulting from the problem Hamiltonian.\"\"\"\n",
    "    return np.array(\n",
    "        [\n",
    "            [np.cos(alpha), -np.sin(alpha), 0],\n",
    "            [np.sin(alpha), np.cos(alpha), 0],\n",
    "            [0, 0, 1],\n",
    "        ]\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "391beb94-5f8b-49f5-8a5e-8df3b4d5294c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def V_D(alpha):\n",
    "    \"\"\"Returns the rotation matrix resulting from the driver Hamiltonian.\"\"\"\n",
    "    return np.array(\n",
    "        [\n",
    "            [1, 0, 0],\n",
    "            [0, np.cos(alpha), -np.sin(alpha)],\n",
    "            [0, np.sin(alpha), np.cos(alpha)],\n",
    "        ]\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7475058d-9209-419e-9b62-12c96d8ad2f7",
   "metadata": {},
   "source": [
    "这个函数是用来计算在演化过程中的所有量子的状态"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "16c608e7-c129-4f13-a901-749c9ec0369b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def evolve(S, h, J, beta, gamma):\n",
    "    \"\"\"\n",
    "    Returns the full history of the vector of spin vectors.\n",
    "\n",
    "    Args：\n",
    "        S (np.ndarray), the initial vector of all spin vectors.\n",
    "        h (list), the vector of local magnetic fields of the problem Hamiltonian.\n",
    "        J (np.ndarray), The coupling matrix of the problem Hamiltonian.\n",
    "        beta (List), the vector of QAOA driver parameters.\n",
    "        gamma (List), the vector of QAOA problem parameters.\n",
    "\n",
    "    Returns：\n",
    "        the full history of the vector of all spin vectors after a full mean-field AOA evolution.\n",
    "    \"\"\"\n",
    "    for k in range(len(beta)):\n",
    "        current_S = S[k, :, :]\n",
    "        m = magnetization(current_S, h, J, len(h))\n",
    "        next_S = []\n",
    "        for i in range(len(h)):\n",
    "            theta = -2 * gamma[k] * m[i]\n",
    "            phi = -2 * beta[k]\n",
    "            v_P = V_P(theta)\n",
    "            v_D = V_D(phi)\n",
    "            new_spin = np.dot(v_D, np.dot(v_P, S[k, i, :]))\n",
    "            next_S.append(new_spin.tolist())\n",
    "        S[k + 1, :, :] = next_S\n",
    "    return S"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf6421c3-12dd-45c3-9547-057cddc7ac57",
   "metadata": {},
   "source": [
    "这个函数是用来计算问题的能量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3a5d877d-c8e1-4718-903a-96f170c2fad5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def expectation(S, h, J):\n",
    "    \"\"\"\n",
    "    Args:\n",
    "        S (np.ndarray), a vector of all spin vectors.\n",
    "        h (list), the vector of local magnetic fields of the problem Hamiltonian.\n",
    "        J (np.ndarray), The coupling matrix of the problem Hamiltonian.\n",
    "\n",
    "    return:\n",
    "        result (float), The energy expectation value corresponding to the supplied spin configuration.\n",
    "    \"\"\"\n",
    "    S_z = np.array([S[i][2] for i in range(S.shape[0])])\n",
    "    result = -(\n",
    "        (np.transpose(h) + 0.5 * np.transpose(S_z) @ J[: len(S), : len(S)]) @ S_z\n",
    "    )\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "04feaed7-c499-44e5-b464-2922cb491e1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_field_solution(S):\n",
    "    \"\"\"Return the solution bitstring computed from a given vector of spin vectors\"\"\"\n",
    "    S_z = [s[2] for s in S]\n",
    "    return np.sign(S_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfe6966e-198b-4d9a-b39d-c2cfe614875a",
   "metadata": {},
   "source": [
    "这个函数是用来计算问题哈密顿量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "61fd24d9-3b03-4726-a99a-16a6e4af0048",
   "metadata": {},
   "outputs": [],
   "source": [
    "def problem_hamiltonian(problem):\n",
    "    \"\"\"Returns the problem Hamiltonian corresponding to 'problem'.\"\"\"\n",
    "    num_qubits = problem.num_qubits\n",
    "    local_fields = problem.local_fields\n",
    "    couplings = problem.couplings\n",
    "    H = QubitOperator()\n",
    "    for i in range(num_qubits):\n",
    "        H += QubitOperator(f\"Z{i}\", coefficient=local_fields[i])\n",
    "    for i in range(num_qubits):\n",
    "        for j in range(i + 1, num_qubits):\n",
    "            H += QubitOperator(f\"Z{i} Z{j}\", coefficient=couplings[i][j])\n",
    "    return H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45204ada-d15e-4ed3-8f05-98170b6a1458",
   "metadata": {},
   "source": [
    "这个函数是用来计算驱动哈密顿量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "80981fe1-7678-4a9b-a461-f36e7d6bb4d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def driver_hamiltonian(problem):\n",
    "    \"\"\"Returns the driver Hamiltonian corresponding to 'problem'.\"\"\"\n",
    "    num_qubits = problem.num_qubits\n",
    "    H = QubitOperator()\n",
    "    for i in range(num_qubits):\n",
    "        H += QubitOperator(f\"X{i}\", coefficient=1)\n",
    "    return H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc3e26f7-7fdd-4068-9b29-eb54e3fac5c2",
   "metadata": {},
   "source": [
    "<span style=\"font-size:17px;\">4 输出结果</span>\n",
    "\n",
    "J矩阵是原论文中的标记为\"easy\"的例子，来实现第一幅图的复现\n",
    "\n",
    "步骤是：\n",
    "\n",
    "1.初始化全体旋转量子的坐标为[1,0,0]\n",
    "\n",
    "2.进行演化迭代，得到每一次迭代返回的量子状态\n",
    "\n",
    "3.计算问题层能量(包括局域场和耦合矩阵)\n",
    "\n",
    "4.计算驱动层的能量(只关于量子的x分量和)\n",
    "\n",
    "5.计算总能量E(t)=-2*(γ(t)*问题层能量+β(t)*驱动层能量) 其中γ(t)是从0到0.5的线性函数  β(t)是从0.5到0的线性函数\n",
    "\n",
    "横坐标是归一化的时间，纵坐标是E-E0，QAOA算法的20条蓝线输出的是QAOA算法前第2-21个能量E-E0，红线输出的是Mean-field AOA算法的能量E-E0，其中E0始终为QAOA的第一能量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ce239652-ecd7-4c83-b834-81fc07a124dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fig7_1(t=0.5):\n",
    "    \"\"\"Use the J in the 'easy' instance\"\"\"\n",
    "    num_layers = 2000\n",
    "    J = np.array(\n",
    "        [\n",
    "            [\n",
    "                0.0,\n",
    "                0.42132292,\n",
    "                -0.25571582,\n",
    "                -0.16267926,\n",
    "                -0.77219702,\n",
    "                0.06310006,\n",
    "                0.03975859,\n",
    "                0.14472543,\n",
    "                0.52842128,\n",
    "                0.2399146,\n",
    "                -0.5426827,\n",
    "            ],\n",
    "            [\n",
    "                0.42132292,\n",
    "                0.0,\n",
    "                -0.04705895,\n",
    "                -0.28587394,\n",
    "                0.50008779,\n",
    "                0.19371417,\n",
    "                0.43070882,\n",
    "                0.21708425,\n",
    "                -0.09432905,\n",
    "                -0.07647162,\n",
    "                0.0410539,\n",
    "            ],\n",
    "            [\n",
    "                -0.25571582,\n",
    "                -0.04705895,\n",
    "                0.0,\n",
    "                0.04119337,\n",
    "                -0.03143317,\n",
    "                -0.69452537,\n",
    "                -0.20650927,\n",
    "                0.03894435,\n",
    "                -0.03114271,\n",
    "                0.31022423,\n",
    "                -0.25880392,\n",
    "            ],\n",
    "            [\n",
    "                -0.16267926,\n",
    "                -0.28587394,\n",
    "                0.04119337,\n",
    "                0.0,\n",
    "                0.46391918,\n",
    "                -0.46363437,\n",
    "                -0.02663793,\n",
    "                0.05954145,\n",
    "                -0.35251299,\n",
    "                -0.13726513,\n",
    "                -0.07634969,\n",
    "            ],\n",
    "            [\n",
    "                -0.77219702,\n",
    "                0.50008779,\n",
    "                -0.03143317,\n",
    "                0.46391918,\n",
    "                0.0,\n",
    "                0.07847923,\n",
    "                0.18372387,\n",
    "                -0.47860755,\n",
    "                -0.24819476,\n",
    "                0.27132491,\n",
    "                -0.22750873,\n",
    "            ],\n",
    "            [\n",
    "                0.06310006,\n",
    "                0.19371417,\n",
    "                -0.69452537,\n",
    "                -0.46363437,\n",
    "                0.07847923,\n",
    "                0.0,\n",
    "                0.01450035,\n",
    "                0.29893034,\n",
    "                0.12659777,\n",
    "                -0.20628129,\n",
    "                -0.17935806,\n",
    "            ],\n",
    "            [\n",
    "                0.03975859,\n",
    "                0.43070882,\n",
    "                -0.20650927,\n",
    "                -0.02663793,\n",
    "                0.18372387,\n",
    "                0.01450035,\n",
    "                0.0,\n",
    "                0.51434209,\n",
    "                -0.13210798,\n",
    "                0.28592817,\n",
    "                0.1470975,\n",
    "            ],\n",
    "            [\n",
    "                0.14472543,\n",
    "                0.21708425,\n",
    "                0.03894435,\n",
    "                0.05954145,\n",
    "                -0.47860755,\n",
    "                0.29893034,\n",
    "                0.51434209,\n",
    "                0.0,\n",
    "                -0.250012,\n",
    "                0.0366756,\n",
    "                -0.59322573,\n",
    "            ],\n",
    "            [\n",
    "                0.52842128,\n",
    "                -0.09432905,\n",
    "                -0.03114271,\n",
    "                -0.35251299,\n",
    "                -0.24819476,\n",
    "                0.12659777,\n",
    "                -0.13210798,\n",
    "                -0.250012,\n",
    "                0.0,\n",
    "                0.26835453,\n",
    "                0.00088238,\n",
    "            ],\n",
    "            [\n",
    "                0.2399146,\n",
    "                -0.07647162,\n",
    "                0.31022423,\n",
    "                -0.13726513,\n",
    "                0.27132491,\n",
    "                -0.20628129,\n",
    "                0.28592817,\n",
    "                0.0366756,\n",
    "                0.26835453,\n",
    "                0.0,\n",
    "                -0.48081058,\n",
    "            ],\n",
    "            [\n",
    "                -0.5426827,\n",
    "                0.0410539,\n",
    "                -0.25880392,\n",
    "                -0.07634969,\n",
    "                -0.22750873,\n",
    "                -0.17935806,\n",
    "                0.1470975,\n",
    "                -0.59322573,\n",
    "                0.00088238,\n",
    "                -0.48081058,\n",
    "                0.0,\n",
    "            ],\n",
    "        ]\n",
    "    )\n",
    "    problem = Problem(num_layers=num_layers, couplings=J)\n",
    "    local_fields = problem.local_fields\n",
    "    couplings = J\n",
    "    num_qubits = problem.num_qubits\n",
    "    beta = [t * (1 - i / num_layers) for i in range(num_layers)]\n",
    "    gamma = [t * (i + 1) / num_layers for i in range(num_layers)]\n",
    "    S = np.zeros((num_layers + 1, num_qubits, 3), dtype=np.float64)\n",
    "    S[0, :, :] = np.array([[1, 0, 0] for _ in range(num_qubits)])\n",
    "    S = evolve(S, local_fields, couplings, beta, gamma)\n",
    "    Ep = [expectation(S[i], local_fields, J) for i in range(num_layers)]\n",
    "    Ed = [-sum([S[i, j, 0] for j in range(num_qubits)]) for i in range(num_layers)]\n",
    "    hamiltonian = [\n",
    "        Hamiltonian(\n",
    "            -2 * gamma[i] * problem_hamiltonian(problem)\n",
    "            - 2 * beta[i] * driver_hamiltonian(problem)\n",
    "        )\n",
    "        for i in range(len(beta))\n",
    "    ]\n",
    "    eig_all = np.zeros((num_layers, 21))\n",
    "    for i in range(len(gamma)):\n",
    "        eigenvalues, _ = eigsh(hamiltonian[i].sparse_matrix, k=21, which=\"SA\")\n",
    "        E0 = eigenvalues[0]\n",
    "        eig_all[i, 0] = E0\n",
    "        for j in range(20):\n",
    "            eig_all[i, j + 1] = eigenvalues[j + 1] - E0\n",
    "    E0 = eig_all[:, 0]\n",
    "    E_meanfield = [2 * gamma[i] * Ep[i] + 2 * beta[i] * Ed[i] for i in range(len(Ep))]\n",
    "    draw_graph1(eig_all, E_meanfield - E0)\n",
    "\n",
    "\n",
    "# t = 0.5\n",
    "# fig7_1(t)\n",
    "# remove the annotations to present FIG 7.1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "488ae4e1-a8fb-4843-9c2a-b2a67cb0ae29",
   "metadata": {},
   "source": [
    "J矩阵是原论文中的标记为\"hard\"的例子，来实现第一幅图的复现\n",
    "\n",
    "步骤：\n",
    "\n",
    "1.初始化全体旋转量子的坐标为[1,0,0]\n",
    "\n",
    "2.进行演化迭代，得到每一次迭代返回的量子状态\n",
    "\n",
    "3.计算问题层能量(包括局域场和耦合矩阵)\n",
    "\n",
    "4.计算驱动层的能量(只关于量子的x分量和)\n",
    "\n",
    "5.计算总能量E(t)=-2*(γ(t)*问题层能量+β(t)*驱动层能量) 其中γ(t)是从0到0.5的线性函数  β(t)是从0.5到0的线性函数\n",
    "\n",
    "横坐标是归一化的时间，纵坐标是E-E0，QAOA算法的20条蓝线输出的是QAOA算法前第2-21个能量E-E0 红线输出的是Mean-field AOA算法的能量E-E0，其中E0始终为QAOA的第一能量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6a23aebc-3962-4387-b205-e307c6908085",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fig7_2(t=0.5):\n",
    "    \"\"\"Use the J in the 'hard' instance\"\"\"\n",
    "    num_layers = 5000\n",
    "    J = np.array(\n",
    "        [\n",
    "            [\n",
    "                0.0,\n",
    "                0.20437418,\n",
    "                0.18668785,\n",
    "                -0.42714492,\n",
    "                -0.49929211,\n",
    "                -0.05025879,\n",
    "                -0.08892105,\n",
    "                -0.39034422,\n",
    "                0.45183176,\n",
    "                0.13354867,\n",
    "                0.15657114,\n",
    "            ],\n",
    "            [\n",
    "                0.20437418,\n",
    "                0.0,\n",
    "                0.47448795,\n",
    "                -0.11285034,\n",
    "                -0.19304397,\n",
    "                -0.54407782,\n",
    "                -0.10369845,\n",
    "                -0.35005243,\n",
    "                0.56644052,\n",
    "                0.07814467,\n",
    "                -0.23853172,\n",
    "            ],\n",
    "            [\n",
    "                0.18668785,\n",
    "                0.47448795,\n",
    "                0.0,\n",
    "                0.09938403,\n",
    "                -0.66523568,\n",
    "                -0.19274242,\n",
    "                -0.08422971,\n",
    "                0.16213894,\n",
    "                -0.03594969,\n",
    "                0.05643456,\n",
    "                0.03460391,\n",
    "            ],\n",
    "            [\n",
    "                -0.42714492,\n",
    "                -0.11285034,\n",
    "                0.09938403,\n",
    "                0.0,\n",
    "                -0.14651126,\n",
    "                -0.57655351,\n",
    "                -0.06858528,\n",
    "                0.02506287,\n",
    "                0.15807228,\n",
    "                -0.05411843,\n",
    "                0.10554402,\n",
    "            ],\n",
    "            [\n",
    "                -0.49929211,\n",
    "                -0.19304397,\n",
    "                -0.66523568,\n",
    "                -0.14651126,\n",
    "                0.0,\n",
    "                -0.42172937,\n",
    "                -0.07213464,\n",
    "                0.19051599,\n",
    "                0.29572682,\n",
    "                -0.46951685,\n",
    "                -0.01001654,\n",
    "            ],\n",
    "            [\n",
    "                -0.05025879,\n",
    "                -0.54407782,\n",
    "                -0.19274242,\n",
    "                -0.57655351,\n",
    "                -0.42172937,\n",
    "                0.0,\n",
    "                0.22313245,\n",
    "                -0.11621411,\n",
    "                -0.58635984,\n",
    "                -0.24611057,\n",
    "                -0.32501641,\n",
    "            ],\n",
    "            [\n",
    "                -0.08892105,\n",
    "                -0.10369845,\n",
    "                -0.08422971,\n",
    "                -0.06858528,\n",
    "                -0.07213464,\n",
    "                0.22313245,\n",
    "                0.0,\n",
    "                -0.3046252,\n",
    "                -0.0152095,\n",
    "                -0.22280432,\n",
    "                -0.43076371,\n",
    "            ],\n",
    "            [\n",
    "                -0.39034422,\n",
    "                -0.35005243,\n",
    "                0.16213894,\n",
    "                0.02506287,\n",
    "                0.19051599,\n",
    "                -0.11621411,\n",
    "                -0.3046252,\n",
    "                0.0,\n",
    "                -0.64309655,\n",
    "                0.24139295,\n",
    "                -0.11473274,\n",
    "            ],\n",
    "            [\n",
    "                0.45183176,\n",
    "                0.56644052,\n",
    "                -0.03594969,\n",
    "                0.15807228,\n",
    "                0.29572682,\n",
    "                -0.58635984,\n",
    "                -0.0152095,\n",
    "                -0.64309655,\n",
    "                0.0,\n",
    "                -0.09709491,\n",
    "                -0.39243752,\n",
    "            ],\n",
    "            [\n",
    "                0.13354867,\n",
    "                0.07814467,\n",
    "                0.05643456,\n",
    "                -0.05411843,\n",
    "                -0.46951685,\n",
    "                -0.24611057,\n",
    "                -0.22280432,\n",
    "                0.24139295,\n",
    "                -0.09709491,\n",
    "                0.0,\n",
    "                0.43852014,\n",
    "            ],\n",
    "            [\n",
    "                0.15657114,\n",
    "                -0.23853172,\n",
    "                0.03460391,\n",
    "                0.10554402,\n",
    "                -0.01001654,\n",
    "                -0.32501641,\n",
    "                -0.43076371,\n",
    "                -0.11473274,\n",
    "                -0.39243752,\n",
    "                0.43852014,\n",
    "                0.0,\n",
    "            ],\n",
    "        ]\n",
    "    )\n",
    "    problem = Problem(num_layers=num_layers, couplings=J)\n",
    "    local_fields = problem.local_fields\n",
    "    couplings = J\n",
    "    num_qubits = problem.num_qubits\n",
    "    beta = [t * (1 - i / num_layers) for i in range(num_layers)]\n",
    "    gamma = [t * (i + 1) / num_layers for i in range(num_layers)]\n",
    "    S = np.zeros((num_layers + 1, num_qubits, 3), dtype=np.float64)\n",
    "    S[0, :, :] = np.array([[1.0, 0.0, 0.0] for _ in range(num_qubits)])\n",
    "    S = evolve(S, local_fields, couplings, beta, gamma)\n",
    "    Ep = [expectation(S[i], local_fields, J) for i in range(num_layers)]\n",
    "    Ed = [-sum([S[i, j, 0] for j in range(num_qubits)]) for i in range(num_layers)]\n",
    "    hamiltonian = [\n",
    "        Hamiltonian(\n",
    "            -2 * gamma[i] * problem_hamiltonian(problem)\n",
    "            - 2 * beta[i] * driver_hamiltonian(problem)\n",
    "        )\n",
    "        for i in range(len(beta))\n",
    "    ]\n",
    "    eig_all = np.zeros((num_layers, 21))\n",
    "    for i in range(len(gamma)):\n",
    "        eigenvalues, _ = eigsh(hamiltonian[i].sparse_matrix, k=21, which=\"SA\")\n",
    "        E0 = eigenvalues[0]\n",
    "        eig_all[i, 0] = E0\n",
    "        for j in range(20):\n",
    "            eig_all[i, j + 1] = eigenvalues[j + 1] - E0\n",
    "    E0 = eig_all[:, 0]\n",
    "    E_meanfield = [2 * gamma[i] * Ep[i] + 2 * beta[i] * Ed[i] for i in range(len(Ep))]\n",
    "    draw_graph1(eig_all, E_meanfield - E0)\n",
    "\n",
    "\n",
    "# t = 0.5\n",
    "# fig7_2(t)\n",
    "# remove the annotations to present FIG 7.2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae14eb96-c54d-4a3c-9175-9cbd9d046b35",
   "metadata": {},
   "source": [
    "J矩阵是原论文中的标记为\"easy\"的例子(同FIG7.1)，来实现第三幅图的复现\n",
    "\n",
    "步骤：\n",
    "\n",
    "1.初始化全体旋转量子的坐标为[1,0,0]\n",
    "\n",
    "2.进行演化迭代，得到每一次迭代返回的量子状态\n",
    "\n",
    "3.用每一次迭代的量子状态计算每一次迭代的lyapunov矩阵\n",
    "\n",
    "横坐标是归一化的时间s，纵坐标是Lyapunov指数，最终输出了Mean-field AOA算法的前三个Lyapunov指数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "93c69a0a-4495-49ed-a2d7-679039a8b2a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fig7_3(t=0.5):\n",
    "    \"\"\"Use the J in the 'easy' instance\"\"\"\n",
    "    num_layers = 2000\n",
    "    J = np.array(\n",
    "        [\n",
    "            [\n",
    "                0.0,\n",
    "                0.42132292,\n",
    "                -0.25571582,\n",
    "                -0.16267926,\n",
    "                -0.77219702,\n",
    "                0.06310006,\n",
    "                0.03975859,\n",
    "                0.14472543,\n",
    "                0.52842128,\n",
    "                0.2399146,\n",
    "                -0.5426827,\n",
    "            ],\n",
    "            [\n",
    "                0.42132292,\n",
    "                0.0,\n",
    "                -0.04705895,\n",
    "                -0.28587394,\n",
    "                0.50008779,\n",
    "                0.19371417,\n",
    "                0.43070882,\n",
    "                0.21708425,\n",
    "                -0.09432905,\n",
    "                -0.07647162,\n",
    "                0.0410539,\n",
    "            ],\n",
    "            [\n",
    "                -0.25571582,\n",
    "                -0.04705895,\n",
    "                0.0,\n",
    "                0.04119337,\n",
    "                -0.03143317,\n",
    "                -0.69452537,\n",
    "                -0.20650927,\n",
    "                0.03894435,\n",
    "                -0.03114271,\n",
    "                0.31022423,\n",
    "                -0.25880392,\n",
    "            ],\n",
    "            [\n",
    "                -0.16267926,\n",
    "                -0.28587394,\n",
    "                0.04119337,\n",
    "                0.0,\n",
    "                0.46391918,\n",
    "                -0.46363437,\n",
    "                -0.02663793,\n",
    "                0.05954145,\n",
    "                -0.35251299,\n",
    "                -0.13726513,\n",
    "                -0.07634969,\n",
    "            ],\n",
    "            [\n",
    "                -0.77219702,\n",
    "                0.50008779,\n",
    "                -0.03143317,\n",
    "                0.46391918,\n",
    "                0.0,\n",
    "                0.07847923,\n",
    "                0.18372387,\n",
    "                -0.47860755,\n",
    "                -0.24819476,\n",
    "                0.27132491,\n",
    "                -0.22750873,\n",
    "            ],\n",
    "            [\n",
    "                0.06310006,\n",
    "                0.19371417,\n",
    "                -0.69452537,\n",
    "                -0.46363437,\n",
    "                0.07847923,\n",
    "                0.0,\n",
    "                0.01450035,\n",
    "                0.29893034,\n",
    "                0.12659777,\n",
    "                -0.20628129,\n",
    "                -0.17935806,\n",
    "            ],\n",
    "            [\n",
    "                0.03975859,\n",
    "                0.43070882,\n",
    "                -0.20650927,\n",
    "                -0.02663793,\n",
    "                0.18372387,\n",
    "                0.01450035,\n",
    "                0.0,\n",
    "                0.51434209,\n",
    "                -0.13210798,\n",
    "                0.28592817,\n",
    "                0.1470975,\n",
    "            ],\n",
    "            [\n",
    "                0.14472543,\n",
    "                0.21708425,\n",
    "                0.03894435,\n",
    "                0.05954145,\n",
    "                -0.47860755,\n",
    "                0.29893034,\n",
    "                0.51434209,\n",
    "                0.0,\n",
    "                -0.250012,\n",
    "                0.0366756,\n",
    "                -0.59322573,\n",
    "            ],\n",
    "            [\n",
    "                0.52842128,\n",
    "                -0.09432905,\n",
    "                -0.03114271,\n",
    "                -0.35251299,\n",
    "                -0.24819476,\n",
    "                0.12659777,\n",
    "                -0.13210798,\n",
    "                -0.250012,\n",
    "                0.0,\n",
    "                0.26835453,\n",
    "                0.00088238,\n",
    "            ],\n",
    "            [\n",
    "                0.2399146,\n",
    "                -0.07647162,\n",
    "                0.31022423,\n",
    "                -0.13726513,\n",
    "                0.27132491,\n",
    "                -0.20628129,\n",
    "                0.28592817,\n",
    "                0.0366756,\n",
    "                0.26835453,\n",
    "                0.0,\n",
    "                -0.48081058,\n",
    "            ],\n",
    "            [\n",
    "                -0.5426827,\n",
    "                0.0410539,\n",
    "                -0.25880392,\n",
    "                -0.07634969,\n",
    "                -0.22750873,\n",
    "                -0.17935806,\n",
    "                0.1470975,\n",
    "                -0.59322573,\n",
    "                0.00088238,\n",
    "                -0.48081058,\n",
    "                0.0,\n",
    "            ],\n",
    "        ]\n",
    "    )\n",
    "    problem = Problem(num_layers=num_layers, couplings=J)\n",
    "    beta = [t * (1 - i / problem.num_layers) for i in range(problem.num_layers)]\n",
    "    gamma = [t * (i + 1) / problem.num_layers for i in range(problem.num_layers)]\n",
    "    lamda = evolve_fluctuations(problem, t, beta, gamma)\n",
    "    lamda_first = lamda[:, -1]\n",
    "    lamda_second = lamda[:, -2]\n",
    "    lamda_third = lamda[:, -3]\n",
    "    draw_graph(lamda_first, lamda_second, lamda_third)\n",
    "\n",
    "\n",
    "# t = 0.5\n",
    "# fig7_3(t)\n",
    "# remove the annotations to present FIG 7.3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32180d90-c76e-4bfd-b43b-76cb568bcbf4",
   "metadata": {},
   "source": [
    "J矩阵是原论文中的标记为\"hard\"的例子(同FIG7.2)，来实现第四幅图的复现\n",
    "\n",
    "步骤：\n",
    "\n",
    "1.初始化全体旋转量子的坐标为[1,0,0]\n",
    "\n",
    "2.进行演化迭代，得到每一次迭代返回的量子状态\n",
    "\n",
    "3.用每一次迭代的量子状态计算每一次迭代的lyapunov矩阵\n",
    "\n",
    "横坐标是归一化的时间s，纵坐标是Lyapunov指数，最终输出了Mean-field AOA算法的前三个Lyapunov指数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "97f82ec8-0268-4586-bd32-442b1dcac074",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fig7_4(t=0.5):\n",
    "    \"\"\"Use the J in the 'hard' instance\"\"\"\n",
    "    num_layers = 5000\n",
    "    J = np.array(\n",
    "        [\n",
    "            [\n",
    "                0.0,\n",
    "                0.20437418,\n",
    "                0.18668785,\n",
    "                -0.42714492,\n",
    "                -0.49929211,\n",
    "                -0.05025879,\n",
    "                -0.08892105,\n",
    "                -0.39034422,\n",
    "                0.45183176,\n",
    "                0.13354867,\n",
    "                0.15657114,\n",
    "            ],\n",
    "            [\n",
    "                0.20437418,\n",
    "                0.0,\n",
    "                0.47448795,\n",
    "                -0.11285034,\n",
    "                -0.19304397,\n",
    "                -0.54407782,\n",
    "                -0.10369845,\n",
    "                -0.35005243,\n",
    "                0.56644052,\n",
    "                0.07814467,\n",
    "                -0.23853172,\n",
    "            ],\n",
    "            [\n",
    "                0.18668785,\n",
    "                0.47448795,\n",
    "                0.0,\n",
    "                0.09938403,\n",
    "                -0.66523568,\n",
    "                -0.19274242,\n",
    "                -0.08422971,\n",
    "                0.16213894,\n",
    "                -0.03594969,\n",
    "                0.05643456,\n",
    "                0.03460391,\n",
    "            ],\n",
    "            [\n",
    "                -0.42714492,\n",
    "                -0.11285034,\n",
    "                0.09938403,\n",
    "                0.0,\n",
    "                -0.14651126,\n",
    "                -0.57655351,\n",
    "                -0.06858528,\n",
    "                0.02506287,\n",
    "                0.15807228,\n",
    "                -0.05411843,\n",
    "                0.10554402,\n",
    "            ],\n",
    "            [\n",
    "                -0.49929211,\n",
    "                -0.19304397,\n",
    "                -0.66523568,\n",
    "                -0.14651126,\n",
    "                0.0,\n",
    "                -0.42172937,\n",
    "                -0.07213464,\n",
    "                0.19051599,\n",
    "                0.29572682,\n",
    "                -0.46951685,\n",
    "                -0.01001654,\n",
    "            ],\n",
    "            [\n",
    "                -0.05025879,\n",
    "                -0.54407782,\n",
    "                -0.19274242,\n",
    "                -0.57655351,\n",
    "                -0.42172937,\n",
    "                0.0,\n",
    "                0.22313245,\n",
    "                -0.11621411,\n",
    "                -0.58635984,\n",
    "                -0.24611057,\n",
    "                -0.32501641,\n",
    "            ],\n",
    "            [\n",
    "                -0.08892105,\n",
    "                -0.10369845,\n",
    "                -0.08422971,\n",
    "                -0.06858528,\n",
    "                -0.07213464,\n",
    "                0.22313245,\n",
    "                0.0,\n",
    "                -0.3046252,\n",
    "                -0.0152095,\n",
    "                -0.22280432,\n",
    "                -0.43076371,\n",
    "            ],\n",
    "            [\n",
    "                -0.39034422,\n",
    "                -0.35005243,\n",
    "                0.16213894,\n",
    "                0.02506287,\n",
    "                0.19051599,\n",
    "                -0.11621411,\n",
    "                -0.3046252,\n",
    "                0.0,\n",
    "                -0.64309655,\n",
    "                0.24139295,\n",
    "                -0.11473274,\n",
    "            ],\n",
    "            [\n",
    "                0.45183176,\n",
    "                0.56644052,\n",
    "                -0.03594969,\n",
    "                0.15807228,\n",
    "                0.29572682,\n",
    "                -0.58635984,\n",
    "                -0.0152095,\n",
    "                -0.64309655,\n",
    "                0.0,\n",
    "                -0.09709491,\n",
    "                -0.39243752,\n",
    "            ],\n",
    "            [\n",
    "                0.13354867,\n",
    "                0.07814467,\n",
    "                0.05643456,\n",
    "                -0.05411843,\n",
    "                -0.46951685,\n",
    "                -0.24611057,\n",
    "                -0.22280432,\n",
    "                0.24139295,\n",
    "                -0.09709491,\n",
    "                0.0,\n",
    "                0.43852014,\n",
    "            ],\n",
    "            [\n",
    "                0.15657114,\n",
    "                -0.23853172,\n",
    "                0.03460391,\n",
    "                0.10554402,\n",
    "                -0.01001654,\n",
    "                -0.32501641,\n",
    "                -0.43076371,\n",
    "                -0.11473274,\n",
    "                -0.39243752,\n",
    "                0.43852014,\n",
    "                0.0,\n",
    "            ],\n",
    "        ]\n",
    "    )\n",
    "    problem = Problem(num_layers=num_layers, couplings=J)\n",
    "    beta = [t * (1 - i / problem.num_layers) for i in range(problem.num_layers)]\n",
    "    gamma = [t * (i + 1) / problem.num_layers for i in range(problem.num_layers)]\n",
    "    lamda = evolve_fluctuations(problem, t, beta, gamma)\n",
    "    lamda_first = lamda[:, -1]\n",
    "    lamda_second = lamda[:, -2]\n",
    "    lamda_third = lamda[:, -3]\n",
    "    draw_graph(lamda_first, lamda_second, lamda_third)\n",
    "\n",
    "\n",
    "# t = 0.5\n",
    "# fig7_4(t)\n",
    "# remove the annotations to present FIG 7.4"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
